import plotly.graph_objects as go
import numpy as np
from scipy.optimize import root_scalar
from IPython.display import Markdown

# Define constants and known values
M_1 = 0.5
x_1 = 2
c_f = 0.0025
D = 0.1
g = 1.4

In a constant area duct the coefficient of friction plays the same role as the change in area of a converging frictionless duct, meaning if the flow is originally subsonic the Mach number will increase and if it is originally supersonic the Mach number will decrease \ This means that at some point the flow will choke (reach M=1), this is useful because there is an isentropic flow relation between the Mach number and the distance to the choked position: $$ \require{color}\frac{4c_{f}}{D}L = \frac{1 - {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2}}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2} + \frac{\gamma +1}{2\gamma} \ ln(\frac{\frac{\gamma +1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}{1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}) $$ This finds the maximum required distance for the flow to choke, then since we know the distance between $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}$ and $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ (where we want to know the Mach number) we can also find the distance it takes for the flow to choke from $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ \ Then using the distance bwetween $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ and $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}=1$ we can solve the above equation for Mach number to find $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ \ However, the equation above does not have an analytical solution so a process such as curve fitting or root finding must be used, in this notebook I use root finding using scipy

# define a function that calculates how far it takes for the flow to choke given the initial mach number 
def max_L(M):
    L = (D/(4*c_f))*((1-M**2)/(1.4*M**2) + (g+1)/(2*g) * np.log(((M**2)*(g+1)/2)/(1+(M**2)*(g-1)/2)))
    return L

# define the length to where the flow chokes and a position we want to find the Mach number 
L_choke = max_L(M_1)
x_2 = x_1 + L_choke/2

x=np.linspace(0, x_1+L_choke+1, 100)

# plot a duct and show where each Mach value is 
fig = go.Figure()
fig.add_scatter(x=x, y=np.e**(-x) + 1, fill ='tonexty', fillcolor='rgba(133, 163, 201, 0.5)', line_color='blue')
fig.add_scatter(x=x, y=-np.e**(-x) - 1, fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='blue')
fig.add_vline(x=x_1, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}$")
fig.add_vline(x=x_1+L_choke, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}=1$")
fig.add_vline(x=x_2, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{2}$")
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

# define the length from M_2 to where the flow chokes
L_2 = L_choke + x_1 - x_2

# define function to find Mach roots 
def mach_function(M, L):
    return (D/(4*c_f))*((1-M**2)/(1.4*M**2) + (g+1)/(2*g) * np.log(((M**2)*(g+1)/2)/(1+(M**2)*(g-1)/2))) - L

# using scipy find the root where L=L_2
M_2_root = root_scalar(mach_function, args=(L_2), x0=M_1, x1=1)
M_2 = M_2_root.root

Markdown(fr'The Mach number at a distance ${round(x_2 - x_1, 3)}m$ from $M_1$ the Mach value is ${round(M_2, 3)}$')

The Mach number at a distance $5.345m$ from $M_1$ the Mach value is $0.589$